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It was recently demonstrated that time-dependent PDE problems can numerically 
be solved with a fully pseudospectral scheme, i.e. using spectral expansions with 



I respect to both spatial and time directions (Hennig and Ansorg, 2009 15|). This 

f-*) I was done with the example of simple scalar wave equations in Minkowski spacetime. 

CN ' Here we show that the method can be used to study interesting physical problems 

O ' that are described by systems of nonlinear PDEs. To this end we consider two 1+1 

dimensional problems: radial oscillations of spherically symmetric Newtonian stars 



and time evolution of Gowdy spacetimes as particular cosmological models in general 
relativity. 
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I. INTRODUCTION 

Many dynamical physical phenomena can be mathematically modelled in terms of (hy- 



CSl ' perbolic) partial differential equations (PDEs). Hence, it is of fundamental interest to study 



properties of PDEs and, of course, to construct their solutions. In the special case of linear 
. PDEs, quite comprehensive theory is available and it is possible to find exact solutions in 
^ I many cases. On the other hand, there are large gaps in the mathematical theory of nonlinear 
CN ' equations and exact solutions are rare exceptions. However, many (if not most) interesting 
processes in nature can be described by nonlinear PDEs, so that it is desirable to have 
methods for the construction of at least approximate solutions in these important cases. 
^ I One approach could be the application of analytical approximation techniques, and another 
possibility are numerical methods. In this article, we focus on this second alternative. 

There is a large variety of numerical methods available that can be used to solve non- 
linear PDEs. If one is interested in highly accurate solutions (in order to be as close to an 
exact solution as possible), one particularly useful approach is (pseudo-) spectral methods. 
In the case of dynamical (time-dependent) problems, traditionally a combination of spec- 
tral methods (with respect to the spatial directions) and finite difference methods (in time 
direction) has been used. We refer the reader to for a comprehensive overview of spec- 
tral methods in general relativity as a particular and interesting field of application. Until 
recently, however, there was almost no experience with spectral expansions with respect to 
space and time^. This changed when it was shown in [l5^ that a fully pseudospectral scheme 
can indeed be successfully applied to hyperbolic PDEs. With the example of simple 1+1 
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^ To the best of our knowledge there is only one article that studies parabolic equations by means of spectral 
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dimensional scalar wave equations in Minkowski spacetime, it was demonstrated in [15(] that 
highly accurate numerical solutions can be obtained in this way. A first application of this 
method is given in j2|, where the general relativistic field equations have been solved in the 
interior region of axisymmetric and stationary black holes. 

In this paper, we present two new applications of the fully pseudospectral scheme to 1+1 
dimensional dynamical problems in the context of Newtonian and Einsteinian gravity. Our 
intention is to demonstrate, with interesting and nontrivial examples, that this approach is 
a very promising method, which seems to be applicable to a large class of nonlinear PDE 
problems. In particular, we show that stable and highly accurate long-term evolutions are 
possible. 

This paper is organized as follows. We start by summarizing the basic idea of the fully 
pseudospectral scheme in Sec. [Tll In Sec. Illlt we consider radial oscillations of spherically 
symmetric stars within Newton's theory of gravitation as a first application of the method. 
Afterwards, the time evolution of particular general relativistic cosmological models (Gowdy 
spacetimes with spatial topology) is studied as a second application in Sec. HVl Finally, 
we conclude with a discussion of our results in Sec. El 



II. DESCRIPTION OF THE NUMERICAL METHOD 

In the following we summarize the ideas of the fully pseudospectral scheme. For more 



details we refer to [15| . The underlying numerical method is a generalization of the method 
by Ansorg et al. for solving elliptic PDEs fi\ to the case of hyperbolic equations. 

Our goal is to construct approximations to functions /i, /2, • • • , /m (depending on a spatial 
coordinate, say x, and on time t) that solve a system of m coupled, nonlinear PDEs subject to 
some initial conditions. Depending on the order of the equations, the initial conditions could 
be prescribed function values and derivatives up to some order at some given instant of time 
t = tmin- The numerical solution should then be calculated in a time interval [tmin, ^max], 
i.e. we assume that we are interested in a compact time domain. Moreover, we assume 
that also x takes on values in a bounded domain only^. In addition to initial conditions, one 
generally expects to have some boundary conditions at the boundaries of the spatial domain. 
These conditions might either come from some physical assumptions or from a particular 
degeneracy of the equations at the boundaries so that regular solutions are only possible if 
appropriate boundary conditions are satisfied. 

Our pseudospectral scheme starts by mapping the physical domain of the PDE problem to 
one or several unit squares by means of a coordinate transformation to spectral coordinates 
0", r G [0, 1] X [0, 1] in each domain. In some cases this transformation is a simple rescaling 
of the physical variables x and t with constant factors. However, in other cases the original 
domain might not even be known from the beginning, i.e. one may also be interested in 
free boundary value problems. (An example is provided by the stellar pulsations to be 
considered in Sec. Illlt where the interior of the star is a sphere of radius R{t). However, 
the function R = R{t) can only be determined together with the unknown fluid variables 



^ In the following examples, this will automatically be the case as a consequence of the physical setting 
(compact interior region of an oscillating star/compact cosmological model). However, one could also start 
from a problem formulated on an unbounded spatial (or/and time) domain and perform a compactification. 



An example for this precedure is the conformal compactification of Minkowski spacetime discussed in 15| . 
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while the hydrodynamic equations are solved and is, therefore, not known a priori.) As a 
consequence, the coordinate transformations may involve unknown functions which have to 
be determined from some additional conditions together with the other unknown functions 
/i, z = 1, . . . ,m. 

In a next step, we need to ensure that the initial conditions at r = (corresponding 
to t = tmin) are satisfied. For that purpose, we write the functions /j(cr, r) in terms of 
functions r) in such a way that the initial conditions are satisfied automatically. As 
an example, consider an initial value problem where the function values at r = are given, 
/j(o", 0) = gi{<j). Then introduce new unknowns /j via 

fi{a, t) = gi{a) + Tfi{(T, r). (1) 

If also derivatives of fi need to be prescribed at the initial surface, then one can choose a 
similar expansion including higher order terms in r, see (l5| . 

Once the mathematical problem is reduced to a PDE problem on unit squares and we 
have taken care of the initial conditions, the main idea of a pseudospectral method is to ap- 
proximate the unknown functions fi in terms of finite linear combinations of some set of basis 
functions. In our scheme we choose Chebyshev polynomials and use an approximation of 
the form 

Ua, r)^Y.H C'^jkTji2a - l)Tfc(2r - 1), z = 1, . . . , m, (2) 

j=0 k=0 

for given spectral resolutions (number of polynomials) = N^- + 1, n^- = Nt- + 1. Note that 
it is a simple task to find approximations for the spatial and time derivatives of fi once we 
have the approximation ([2]). 

Now we have to determine the Chebyshev coefficients Cijk- To this end, we introduce a 
suitable set of collocation points. Here, we choose Gauss-Lobatto nodes^, which are defined 
by 

""^^ (^) ' ^'"'''''(^)' •^■ = 0'---'^- k = 0,...,N^. (3) 

Then we impose the condition that our approximations to the functions fi, together with 
the corresponding approximations to their derivatives, satisfy the PDEs (plus the boundary 
conditions, if applicable) exactly at these points. This provides us with a large algebraic 
system of equations for the Chebyshev coefficients Cijk or, equivalently, for the function 
values fi{(Tj, Tk) — and this is the formulation that we actually use. 

Finally, starting from some initial guess, we solve this system iteratively with the Newton- 
Raphson method. This part of the algorithm makes the method computationally expensive, 
since inversion of relatively large matrices is required. However, in the present case of 1+1 
dimensional problems and for a moderate number of unknown functions fi and gridpoints 
{ncryfir ~ 30) this is no problem and the method will be reasonably fast. 



^ It is well- known that a nearly optimal polynomial approximation is achieved for Gauss-Chebyshev collo- 
cation points (Chebyshev roots), see, e.g., [sj. The Gauss-Lobatto points used here (Chebyshev extrema 
plus boundary points) lead to slightly less accurate approximations. However, they have the advantage 
to include gridpoints at the boundaries of our unit squares, which is useful for implementing boundary 
conditions. 
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There are quite different ways of obtaining an initial guess, for which the Newton- Raphson 
method will converge. In some cases it is possible to start from the trivial solution. Some- 
times, one can introduce a parameter, say a, into the problem, such that the solution is 
known for a = and the relevant solution is obtained for a = 1. Then, starting from the 
case with a = 0, one can go to a = 1 in a few steps and always use the solution from the 
previous step as an initial guess for the next one. This procedure is used for the calculation 
of oscillating stars in Sec. IIIIl where we can start from initial data corresponding to a star 
in equilibrium. The finally interesting star (with initial data that significantly deviate from 
equilibrium data) is then reached in about one to three steps. It is also possible to start 
from an analytically obtained approximate solution, for example by using Taylor expansions 
about the initial hypersurface. This method is used in Sec. [IV] for the calculation of cos- 
mological models. Finally, one could also solve the problem with, e.g., a finite difference 
method first and interpolate the resulting numerical values onto the spectral grid to obtain 
an initial guess. 

The described pseudospectral method is a highly implicit numerical scheme^, in which 
the solution at rio- points in spatial direction times Ur points in time direction is obtained 
simultaneously. This is in contrast to the approach in most other numerical PDE methods, 
where the solution is usually calculated at one time slice after the other. However, this 
implicit character plays a crucial role for obtaining highly accurate solutions and is, therefore, 
an essential part of the fully pseudospectral scheme. 

III. NEWTONIAN STELLAR OSCILLATIONS 

We come now to the main part of this paper and study applications of the fully pseu- 
dospectral scheme. As a first example, we consider nonlinear radial stellar oscillations within 
Newton's theory of gravitation. 

A. The physical model 

Variable stars are a class of astrophysical objects whose luminosity (or, more precisely, 
apparent magnitude as seen from Earth) changes over time. This fascinating phenomenon 
has already been intensively studied for a long time, both observationally and theoretically. 
Among the results is the famous period-luminosity relation [13] for classical Cepheids. The 
idea that, for some types of variable stars, the luminosity changes are caused by stellar pul- 
sations dates back almost 100 years to work by Plummer [22t] and Shapley j23[^. Nowadays, 
stellar pulsations are used as an accurate tool to study many aspects of stellar physics. In 
particular, the numerical simulation of stellar oscillations is an important element of the 
investigation of properties of variable stars. In recent times, oscillations have also attracted 
much attention in the context of compact astrophysical objects such as neutron stars, in 
particular in the field of gravitational wave physics. For an overview over the Newtonian 
description of stellar pulsations we refer to 0]. An interesting recent numerical study (based 
on a finite difference method) of general relativistic radial pulsations can be found in 



^ With the term "highly imphcit" we want to emphasize that the numerical values at all (instead of, perhaps, 

only two or three) time slices are coupled in our method. 
^ Interesting details on the history of the theory of stellar pulsations can be found in 12[ . 
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Since our main interest is in the numerical method, we neglect the effects of general 
relativity and study here the simpler equations of Newtonian physics. To this end, we model 
an oscillating star as a fluid ball of variable radius R = R{t), described by the following 
equations. 



AU = AnGp 
P^=P[|r + (v-V)v] =-Vp 
+ V ■ (pv) = 



dp 
at 



(Poisson equation), 
pVU (Euler equation), 

(Continuity equation) 



(4) 
(5) 
(6) 



where U, p, p, v and G denote the gravitational potential, the mass density, the pressure, 
the velocity field and the gravitational constant, respectively. In addition, we choose a 
polytropic equation of state 



p = Kp'^, 7 = 1 + 



1 



n 



(7) 



with polytropic index n = 1 (7 = 2), where K is the polytropic constant. 

In order to obtain a 1 + 1 dimensional problem, we assume spherical symmetry, where all 
functions depend only on a radial coordinate r and on time t. This results in the simplified 
system of equations 

2 

Urr + -Ur = ^TlGp, (8) 

r ' 

-P,r - pU,r, (9) 

+ -pv = 0, (10) 



p(w t + VV^r) 

P,t + {pv), 



where v is the radial velocity and the subscripts denote derivatives with respect to the 
indicated variables. Finally, after eliminating the pressure with ([7]), we arrive at a system of 
three coupled PDEs for the three functions [/, p and v. Since the equations (|9]) and (fTOj) are 
first order in time, we can prescribe the initial values for v and p. But of course we cannot 
choose an initial potential U , since U is already determined by p from ([8]). 

Regular functions U and p inside the star have the property of being even in r, whereas 
V is an odd function of r. Hence, if we replace v with 



w(r, t) :- 



v{r, t) 



(11) 



then we obtain a set of three even functions, U, p, w, which can be considered as functions 
of and t. This will be used when we introduce adapted spectral coordinates in the next 
subsection. 

For an initial configuration at time t = we choose a slightly disturbed equilibrium 
star. In equilibrium (no time-dependence), the above equations for a polytropic star can be 
reduced to the Lane-Emden equation For n = 1, the exact solution is known and we 
obtain the corresponding equilibrium solution 



Peqi 



Pc smc I 7r 



f/eq(r) = ~2Kp, 



1 + sine I TT 



Weq(r) = 0, 



(12) 



where pc is the central mass density, '^o = y |§ is the equilibrium radius and sincx = 
One possibility of disturbing this equilibrium is to consider a non-vanishing initial velocity 
field, i.e. we can start from the initial data 



p(r, 0) = peq(r), U{r, 0) = U,^{r), w{r, 0) = Wo{r) ^ 



(13) 
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for some function Wq = Wo{r). 

Finally, we note that the exterior vacuum region around an oscillating star is trivially 
described by the potential U = —GM/r, where M is the (constant) mass of the star. Hence, 
we can restrict our attention to the interior region. The boundary between both regions 
(i.e. the rim of the star) is characterized by vanishing pressure, p = 0, and, therefore, by 
p = for our equation of state (JTj). 

B. Numerical details 

In order to study the time evolution of a star over a longer period of time (covering several 
oscillation periods), we divide the time interval [0,tmax] into subintervals 

[to,ti],[ti,t2],...,[tjv-i,tw] with ti=tAt, At:=^. (14) 

In this way we can ensure that the unknown functions do not oscillate too strongly in each 
domain. Otherwise, Chebyshev polynomials of much higher degree would be required for 
accurate approximations. 

For the first time domain, initial conditions are given at t = as described in the previous 
subsection (perturbations of the Lane-Emden equilibrium solution, cf. Eq. (fT3l) ). When the 
functions in the first domain have been computed, new initial data are read off for the 
second domain, and so on. Hence, the time domains are evaluated step by step, one after 
the other. (In a future version of the code, a simultaneous treatment of several domains 
might be considered to further increase the numerical accuracy.) 

Spectral coordinates (a, r) G [0, 1] x [0, 1] in each domain i = 0, . . . , N — 1 are now 
introduced via 

r = y/aR{t), t = ti + TAt, (15) 

where the square root takes into account that all the unknowns are regular functions of 
and t as discussed above. This transformation contains the unknown radius R = R{t), which 
can be determined from the surface condition p{r = R{t),t) = 0. 

In terms of these new coordinates, the fluid equations fl8l)- f[TQl) take the following form, 

2aU^aa + 3f/,, = 2'kGR^p, (16) 
f^w,r + 2a (^Rw - w,^ + Rw^ + j^ {2Kp^^ + f/,.) = 0, (17) 

^p,, + 2a(w- p,^ + p{2<yw,, + 3^) = 0. (18) 

Note that these equations contain only spatial derivatives of U , but not U itself. Hence, U 
is determined only up to an additive function of r (or, equivalently, of t). This is physically 
reasonable, since only the potential difference (work) or the gradient (force) are measurable 
quantities. Here, we fix this degree of freedom in the usual way by choosing t/ — )■ at 
infinity. As a consequence, we have the typical 1/r-potential of spherical symmetry in the 
exterior region. Since every 1/r-potential satisfies rU^r + t/ = 0, this condition must also 
hold for the interior potential at the boundary r = i?, as a consequence of a continuous and 
different iable transition at this boundary. Consequently, we impose this equation instead of 
Eq. f lT6|) at the rim of the star. In terms of the spectral coordinates, this boundary condition 
becomes 

2aU,„ + t/ = 0. (19) 
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(a) (b) 




t [ms] t [ms] 



FIG. 1: Two examples of spectral evolutions (solid curves) with different initial data. The dotted 
curve in the left panel shows for comparison the result as obtained with a finite difference code. 

(a) : central density for v{r,0) = 300^ x sin(7rr/ro), 40 time domains with resolution 16 x 18 

(b) : radius for v{r,0) = 4116^ x r/rg, 150 time domains with resolution 18 x 20. 



Alternatively, one could also prescribe the values of U at some timelike curve, for example on 
the symmetry axis, e.g. by setting U to zero there. But this would introduce an "artificial" 
time dependence of the exterior potential. 

The initial conditions for p and w are implemented by replacing these functions in terms 
of functions p and w as illustrated in Eq. ([1]). Therefore, we finally arrive at the set of 
unknown quantities U{cr,T), p{a,T),w{a,T), R{t), which can be calculated by solving the 
three fluid equations together with the boundary condition for U plus the extra condition 
of vanishing surface pressure. 



C. Numerical results 

In the following, we consider examples of Newtonian stars as simple models for neutron 
stars, i.e. we choose parameters leading to stars with radii about 10 km and masses about 
1.4 Mq. To this end, we start from the equilibrium density p(r, 0) = Pcq{r) with initial 
central density pc = 1.9891 x 10^^^. As discussed above, the equilibrium is then distorted 
by adding an initial velocity field w{r,0). Furthermore, we choose a polytropic constant of 
K = 4.5 X 10~^ so that a star in equilibrium would have a radius of tq = 10.29 km. Only 

in the first calculation (Fig. [T^) we use a slightly different value, namely K = 4.25 x 10"'^ i^^z- 
In a first numerical example, a star is evolved for 2.5 ms (corresponding to 10 fundamental 
oscillations), see Fig. [1^, and the result is compared to data obtained with a finite difference 
code^. The diagram shows that both curves agree well, which is a first confirmation that 
the pseudospectral method produces correct results. 

As second example, we evolve a star with different initial conditions for 12 ms (48 fun- 
damental oscillations). The result is shown in Fig. [lb. Besides the fundamental oscillation 



^ I am grateful to Ernazar Abdikamalov for providing me with the finite difference result. 
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(with frequency^ /o = 4.009596 kHz, corresponding to a period of 0.249402 ms) one can 
see an amplitude modulation with lower frequency. This happens since the frequencies of 
different oscillation modes are not precisely multiples of the fundamental frequency. As a 
consequence, there is a beat with frequency /i — 2/o = 0.25 Hz (period of 3.99 ms), where 
/i = 8.269880 kHz is the frequency of the first overtone. In addition, the amplitude in- 
creases due to another beat with frequency /2 — 3/o = 0.035 kHz (period of 28.2 ms), where 
/a = 12.064215 kHz. 



As illustrated with the numerical examples in the previous subsection (in particular, by 
comparing with the finite difference result), a fully pseudospectral scheme is able to produce 
correct time evolutions of oscillating stars. However, the main reason for using a spectral 
method is, of course, that one wants to produce very accurate results. Therefore, it is desir- 
able to test convergence properties of the code and to measure the accuracy quantitatively. 

To this end, we look at two different quantities and their accuracy. In a first step, we 
may calculate the mass M of the star as a function of time. 



The integral in the latter formula can be approximated by using the Chebyshev represen- 
tation of p(r, t) with the numerically calculated coefficients and integrating the Chebyshev 
polynomials exactly. As a consequence of the continuity equation ([6]) [or ffTOj) ]. M must, 
of course, be independent of t. On the other hand, numerical errors will introduce some 
artificial time dependence of M . This effect can be used as a measure for the numerical 
accuracy. 

The relative mass change during ten oscillations of a star is plotted in Fig. [2^. As one 
might have expected, the average error is growing with time. Finally, the error reaches values 
between 10~^^ and 10~^° at the end of the simulation, which shows that the numerical result 
is indeed very accurate. 

Next, we investigate how the error depends on the spectral resolution no-, Ur and on the 
number of time domains. For that purpose, we look at the numerical values for the central 
mass density at t = tmax as obtained for different resolutions. These values may then be 
compared to the most accurate value (corresponding to the highest resolution), which can be 
considered to be more or less exact. In this way, we obtain a measure for the relative error 
in the central density at the end of the simulations (where the largest errors are expected). 
The result can be found in Fig. Obviously, the error trends downwards for increasing 
spectral resolution, if the number of time domains is fixed. On the other hand, the error 
for given resolution in each domain decreases, if the number of time domains is increased. 
However, once the error has reached values in the order of 10^^^ to 10^^'^ (close to machine 



D. Convergence 



m 




(20) 







^ The numerical values for the frequencies of the oscillation modes have been obtained from a linear per- 
turbation analysis with a pseudospectral code. For details about linear perturbations of polytropic stars, 
see, e.g., or 19|. 
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(a) (b) 




t [ms] spectral resolution ria = Ut in each domain 



FIG. 2: Test of the numerical accuracy: A star with initial velocity profile v{r, 0) = 3087^ x r/rg 
is evolved for 2.5 ms (10 fundamental oscillations). The figure shows: 

(a) The relative mass change \ [m{t) — m(0)]/m(0)| for 25 time domains with a resolution of 25 x 25. 

(b) The relative error in the final central mass density /9(0,tinax) for different numbers of time 
domains and spectral resolutions as estimated by comparing the values of /3(0, tmax) with the value 
obtained for 30 domains with resolution 30 x 30. 



accuracy for our double precision code), it cannot be reduced any further by considering 
more domains or higher resolution. 

This demonstrates that highly accurate solutions can he obtained for moderate values 
of the resolution and number of time domains. Moreover, the error curves in Fig. [2]d are 
approximately linear in a logarithmic plot, i.e. we observe exponential convergence. 



E. Eigenmodes and nonlinear effects 

In order to study effects coming from the nonlinearities in the equations for stellar pulsa- 
tions, we look at the differences between oscillations with small amplitudes and oscillations 
with large amplitudes. 

If we choose initial conditions close to equilibrium (i.e. with a small initial velocity field), 
then the stellar pulsations can be treated as linear perturbations of stars in equilibrium, cf. 
0, EH. As a result, we can expect that the time dependent quantities will undergo harmonic 
oscillations with particular frequencies (corresponding to the eigenmodes of the linearized 
equations). It should also be possible to reproduce this effect with our fully nonlinear code, 
provided the initial perturbations are small enough so that nonlinear terms are irrelevant. 
This can be demonstrated by choosing eigenfunctions of the linear oscillation modes with 
small amplitude as initial velocity field. 

As numerical examples we consider the fundamental mode and the first overtone. The 
resulting radius functions R{t) are shown in Figs. [3^1' andS^. As expected, R{t) seems indeed 
to be a purely harmonic function of time. This can be illustrated more explicitly by looking 
at the corresponding Fourier spectrum as obtained with a discrete Fourier transform, see 
Figs. |3]d and |Dd. Normally, in such a numerical simulation, one might expect to see a large 
peak at the corresponding eigenfrequency together with many superposed peaks (numerical 
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FIG. 3: Fundamental mode with maximal initial velocity v{ro,0) = 20^. The plots show (a) the 
radius R{t) and (b) the Fourier spectrum. 
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FIG. 4: First overtone with maximal initial velocity v{ro,0) = 20— as in Fig. [3l 



noise) due to errors introduced by the numerical method. Here, however, we find nicely 
"clean" spectra, indicating that we have almost pure sinusoidal oscillations. 

Having discussed the linear regime, it is interesting to study larger deviations from equi- 
librium. In that case, nonlinear effects are not negligible anymore. An interesting nonlinear 
effect is the excitation of modes that are not present in the initial data. This can be demon- 
strated by starting from initial data corresponding to the fundamental mode, but with 
increasingly large amplitudes. For small amplitudes, we expect to find again the harmonic 
oscillations with the fundamental frequency. However, for higher amplitudes, we should 
observe that other peaks appear in the spectrum (whose frequencies correspond to higher 
overtones or to linear combinations of several eigenfrequencies) . And indeed, this effect 
shows clearly in the numerical examples, see Fig. [51 
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10 15 20 

frequency [kHz] 

FIG. 5: Starting from data corresponding to the fundamental mode, also other modes are excited 
due to nonlinear effects. The diagrams show the spectra, obtained from Fourier transformation 
of the central mass density p{0,t), for three simulations with different amplitudes of the initial 
velocity field. 

IV. COSMOLOGICAL MODELS 

As a second application of the fully pseudospectral scheme we study a particular cosmo- 
logical model, namely Gowdy spacetimes with spatial topology. 



A. The physical model 



Gowdy spacetimes [13| are relatively simple, but nontrivial inhomogeneous cosmological 
models within Albert Einstein's theory of general relativity. Note that the term "cosmolog- 
ical model" does not imply that we try to find an accurate description of our real universe. 
Instead it refers to a model in mathematical relativity that is assumed to reflect only some 
aspects of the cosmos we live in. A major motivation for the study of these spacetimes 
comes from the desire to understand the mathematical and physical properties of cosmolog- 
ical singularities. 

Mathematically, Gowdy spacetimes are four- dimensional spacetimes characterized by ex- 
istence of two spacelike and commuting Killing vectors (i.e. by existence of two symmetries)^. 
This allows us to introduce adapted coordinates in which everything effectively depends only 
on two (instead of four) coordinates 6 and t. Hence, we again arrive at a 1 + 1 dimensional 
problem. Moreover, we consider vacuum solutions to Einstein's field equations (i.e. no mat- 



® In addition, the precise definition of Gowdy spacetimes also includes that the so-called twist constants 
are zero. 
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ter and no electromagnetic or other fields are present) with vanishing cosmological constant. 
Under these and a few further assumptions^, it turns out that the spatial topology (the 
topology of slices t = constant) must be either T^, S^, S"^ x or that of the lens space 
L{p,q), cf. [H] and [13, H, [2l| . Here, we restrict ourselves to the case of 5*^ topology. 



In order to obtain Gowdy solutions, we have to solve Einstein's field equations. In our 
setting, the essential part of these equations turns out to be equivalent to a single, complex 
equation, the Ernst equation 

^{S) {~£,u - cot tS^t + Em + coteS^e) = -{£,t? + {£,ef (21) 

for the complex Ernst potential £ = £{0,t) =: f{0,t) + ib{6,t). This potential contains 
basically all information about the geometry of the corresponding Gowdy spacetime. 

An important question is whether 5*^ Gowdy solutions that are regular at some time 
t = to will also be regular in the past and future of to- An answer was provided in [g], where 
it was shown that regular data at some to G (0, tt) guarantee regularity in the entire domain 
t G (0, tt), 6* G (0, tt). The boundaries t = and t = tt of this domain are expected to contain 
either curvature singularities or so-called Cauchy horizons^^. Furthermore, the boundaries 
^ = and 6 = n are regular up to coordinate singularities (axes singularities very much like 
the well-known axes singularities of spherical coordinates). 

Here, we consider those Gowdy models with a regular Cauchy horizon at the initial slice 
t = 0. Properties of these spacetimes have been studied in ^ and a family of exact solutions 
has been constructed in The main result is that these Gowdy solutions automatically 
have a second Cauchy horizon at the final slice t = n, unless the initial Ernst potential 
satisfies one of the two conditions 



6(7r,0) -6(0,0) = ±4. (22) 

In this case, the spacetime develops a curvature singularity at t = tt, = (for a sign 
in the latter equation) or at t = vr, ^ = vr (for a '— ' sign), which leads to a divergent Ernst 
potential at the respective points. 

In the following, we intend to solve the Ernst equation (!2T|) numerically, starting from 
initial data at t = 0. Note that, even though (1211) is second order in time, we can only 
prescribe the initial values but not the time derivative of the Ernst potential. The reason is 
that (ETi) degenerates at the special hypersurface t = and leads to the regularity condition 
£_t(^,0) = 0. Moreover, as a consequence of the topology and the fact that t = is 
a Cauchy horizon, the initial potential £{9,0) = So{0) = fo{9) + ibo{9) must be a smooth 
function subject to the following conditions 

/o(0) = /o(7r) = 0, fo{9) > for e G (0, tt), (23) 
/o,e(0) = foA^) = 0, foM^) = /o,ee(7r), (24) 
W(0) = 6o,e(7r) = 0, boMO) = boM^) = 2- (25) 



^ We assume compact, connected, orientable and smooth three manifolds. 

° A Cauchy horizon is a light-hke boundary of the domain in which the Cauchy problem for the Einstein 
equations is valid. The region beyond a Cauchy horizon contains closed time-like geodesies, which implies 
a breakdown of causality and the possibility of time travel backwards in time. 
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B. Numerical details 

According to the above discussion, we expect to find a regular solution to the Ernst 
equation ( !2T1) with initial data at t = in the domain t) G [0, vr] x [0, vr), whereas the slice 
t = IT might contain singularities. Therefore, we intend to construct numerical solutions in 
the domain [0, vr] x [0, tmax] for some tmax € (0, vr). In the singular case, it is then particularly 
interesting to find out how close tmax can be chosen to t = n. 

The next step is to map this domain to unit squares. Here, we choose to represent the 
entire domain by a single square (in contrast to the larger number of time domains as used 
in the example of stellar pulsations). For that purpose, we might consider a simple linear 
rescaling of the coordinates 9 and t. However, it turns out to be useful to choose cos 9 and 
cost as new coordinates first and then to perform an appropriate rescaling, 

1 ^^ 1 — cost 

a = -{1- cos 9), r= , (26) 

because this leads to a particularly nice form of the Ernst equation without trigonometric 
functions. 

It is then straightforward to decompose (1211) into its real and imaginary parts to obtain 
a nonlinear system of two coupled equations for the two unknowns / = and b = '^S. 
Afterwards, the remaining steps of the fully pseudospectral scheme can be carried out with- 
out difficulty as described in Sec. HIl The only point that requires some careful attention is a 
degeneracy of the Ernst equation at the points (0, 0) and (vr, 0), where the reformulation of 
(12T|) in terms of the spectral coordinates a, t turns out to be satisfied identically. However, 
one can derive additional non-trivial conditions at these points from the Ernst equation 



itself. See [15|| for more details about the treatment of such "exceptional points" . 



C. Numerical results 

We consider two numerical examples of solving the Ernst equation (pT!) with the fully 
pseudospectral scheme. 

In the first example, we choose initial data that do not satisfy the characteristic equation 
for singularities (!22|) . Then we know that the solution remains regular in the whole time 
interval [0,7r], and both of the boundaries t = and t = it contain Cauchy horizons. The 
numerical result is illustrated in Fig. |6l As expected, we obtain regular functions 3fJ£ and 
Q£. 

Before we study the accuracy of this result (see Sec. lIVPp . we consider, as a second 
example, initial data subject to (12 2 p with a plus sign on the right hand side. Consequently, 
we expect that the time evolution of these data leads to the formation of a singularity at 
t = IT, 9 = 0. The corresponding numerical solution, shown in Fig. [TJ indeed indicates 
precisely this behaviour. 



D. Convergence 

As in the case of stellar pulsation, a crucial point is the question of convergence and 
numerical accuracy. In order to answer this question in the present case of Gowdy 
spacetimes, we can make use of results from |3j]. In that article, explicit formulae for the 



14 




FIG. 6: Real and imaginary parts of the Ernst potential £ = f + \b obtained from time evolution 
of the initial potential /((9, 0) = sin^O, b{9, 0) = 4(^/3 + 3)V2 + cos6' - 2(^3 + 2) cos d. For these 
initial values we have b{Tr, 0) — 6(0, 0) = 8 — 4-v/3 7^ ±4, which corresponds to a regular case without 
singularities. 




FIG. 7: Numerical values for the real and imaginary parts of the Ernst potential £ = f + ib for 
the initial value problem with f{0,0) = 3sin^^, b{6,0) = cos9 {cos'^6 — cos6 — 3). Here, we have 
b{'ir, 0) — 6(0, 0) =4. As a consequence, the corresponding solution is singular and blows up at 
t = vr, ^ = 0. Interestingly, it is even possible to construct the exact solution to this particular 
initial value problem for the Ernst equation, which is demonstrated in [3]. 



Ernst potential at the boundaries 6 = 0, 6 = ir and t = tt in terms of the initial potential at 
t = have been derived. In particular, for £^ at 6' = we have 

i[&(0, 0)- 2(1 -cost)]g(t,0) + 6^(0,0) 
^ ' ' f(t,0)-i[6(0,0) + 2(1 -cost)] ' ^ ^ 

where S{t,0) means S{6 = t,t = 0). Using this formula, we can, for example, calculate the 
exact value of / = at 9 = 0, t = tmax and compare it with the numerical value. The 
corresponding relative error may then be used as a measure for the numerical accuracy. 

In the following, we apply this method for studying the accuracy of the regular and the 
singular example evolutions presented in the previous subsection. In both cases, we consider 



15 



a number of different spectral resolutions and maximal times tmax- 

The accuracy for the regular example is shown in Fig. [8^. We see that the error decreases 
with increasing resolution until saturation between 10~^^ and 10~^^ is reached. Hence, the 
numerical solution in the regular case is highly accurate. By comparing the error curve for 
^max = 1-5 with the curve for tmax = 3.0, we see that, for larger tmax, a larger resolution 
is required in order to reach the saturation level. But this is clear since evidently more 
Chebyshev polynomials are required if a function is to be accurately approximated on a 
larger domain. 

We can do the same investigation for the singular example studied in the previous sub- 
section. Since the corresponding Ernst potential diverges at t = vr, 6* = 0, we cannot expect 
very accurate solutions if tmax approaches vr. Indeed, this expected drop in the accuracy can 
clearly be seen in Fig. |8]d. For tmax = 1-5, the error decreases to about 10~^^ for n„ = 20 
(and tIt- = 20). In order to reach the same level of accuracy for t^ax = 2.2, we already need 
a resolution of = 30 (and ra^ = 60). Finally, for tmax = 2.6 the error has "only" decreased 
to about 10~^ for the maximal considered resolution of n^^ = 30 (and rij. = 60). This is 
still very accurate compared to errors in typical low-order finite difference algorithms, but 
definitely much below the "usual spectral accuracy" (which is close to machine accuracy). 
Clearly, a huge number of Chebyshev polynomials is needed to approximate functions in the 
vicinity of a singularity. 

A possible approach for increasing the numerical accuracy near singularities could be 
to split the time domain [0, tmax] into several subdomains (as done in the case of stellar 
pulsations). Then, in the first domains a low resolution would be sufficient, whereas 
larger resolutions could be used in the domains close to the singularity. In this way, one 
could "concentrate" higher resolutions where they are needed. Another possibility would 
be to rescale the coordinates in such a way that steep gradients near the singularity are 
reduced. An example for this procedure can be found on p. 119 of jisj. 

Finally, we note that the error curves in the logarithmic plots in Fig. |8] are more or 
less linear, i.e. we again observe exponential convergence (with higher convergence rates for 
regular solutions and lower convergence rates in singular cases). 

V. DISCUSSION 

We have studied the applicability of the fully pseudospectral numerical scheme, intro- 
duced in [l5|, to solving time dependent physical problems. In particular, we have consid- 
ered stellar pulsations of spherically symmetric, polytropic Newtonian stars and the time 
evolution of Gowdy spacetimes with spatial topology. 

The basic idea of the method is to solve initial-boundary value problems for nonlinear 
PDEs with Chebyshev expansions in spatial and time directions. As a consequence, space 
and time are treated relatively equal. Moreover, we obtain a highly implicit method for 
which the numerical solution in entire spacetime domains is obtained simultaneously. 

Some details of the numerical implementation were different in our two examples (free- 
boundary problem and several time domains for stellar oscillations, but fixed boundaries 
and only a single time domain for Gowdy spacetimes). However, the final outcome was the 
same and confirmed the observations from the examples of the simple scalar waves equation 
discussed in [l5||: with a fully pseudospectral scheme, we obtain highly accurate numerical 
solutions with relative error in the order 10~^^ to 10"^'^ (corresponding to up to 13 significant 
digits, which is close to the machine accuracy of our double precision code). This is possible 
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FIG. 8: Test of the numerical accuracy: Two sets of initial data are evolved for different choices 
of the maximal time imax- Afterwards, the numerical value for the real part of the Ernst potential 
at t = tmax, ^ = can be compared with the exact value from the analytical formula (j27p . The 
diagrams show the relative error in the numerical value. 

(a) Initial data as in Fig. El corresponding to a solution with a regular Cauchy horizon at t = vr. 

(b) Initial data as in Fig. [TJ leading to a singularity at t = tt, 9 = 0. 



with moderate spectral resolutions, provided the solutions are sufficiently regular. On the 
other hand, close to singularities, much higher resolutions are required for highly accurate 
solutions; otherwise, the accuracy can drop significantly. However, both in the regular and 
singular cases, the method possesses a geometric convergence rate and the error decreases 
exponentially with the spectral resolution. Moreover, the example of stellar pulsations has 
demonstrated that stable long-term evolutions are possible, namely evolutions over many 
oscillation periods. 

These results should encourage the application of the fully pseudospectral scheme to wide 
classes of problems in physics and other areas. The high accuracy of the method will allow 
for very precise numerical studies, very clean spectra, accurate preservation of conserved 
quantities, etc. 

In the future, it would be desirable to develop the method further in order to allow for a 
treatment of higher dimensional problems (e.g. time evolution of axisymmetric rather than 
spherically symmetric stars). An important step in this direction would be to replace the 
computationally most expensive part of the method, namely the matrix inversion in the 
Newton Raphson scheme, by an iterative method^^. 



The complexity of the matrix inversion via LU-decomposition as used here is proportional to the third 
power of the matrix dimension, and the matrix dimension itself is proportional to x n^.. Moreover, in 
the examples discussed here, about two to four matrix inversions are required in each spectral domain. 
Hence, for ria ^ rir = n, the complexity of the fully pseudospectral scheme is O(n^). As an example, on 
my PC, the inversion of the matrix appearing in the calculation of the Newtonian star takes about 0.18 s 
for n = 15 and about 7.8 s for n = 26. For comparison, a low-order finite difference algorithm would only 
have a complexity of 0{n'^), but, on the other hand, require a huge number of gridpoints to achieve an 
accuracy comparable to "spectral accuracy" . Since a small n is sufficient in the pseudospectral scheme, 
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